Generalized Poland-Scheraga model for supercoiled DNA 
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The Poland-Scheraga (PS) model for the helix-coil transition of DNA considers the statistical 
mechanics of the thermally induced binding of two complementary strands of DNA. In this paper, 
we show how to modify the PS model when a torque is applied to the extremities of DNA: We 
propose a simple model for the energy of twisted DNA and compute the entropy of a loop, subject 
to angular constraints (supercoiling). The denaturation curves are shifted towards lower or higher 
temperatures depending on the sign of the torque, and the UV absorption peaks are softened. The 
properties of supercoiled DNA can be deduced through the use of a numerical Legendre transform. 
In the homogeneous case, we find that for weak supercoiling, the phenomenological quadratic law 
relating the torsional energy to the number of unpaired bases is recovered. 

PACS numbers: 87.14.Gg; 87.15. Cc; 82.39.Pj 

Natural DNA exists as a double helix bound state 0. Upon heating, the two complementary strands may separate. 
This thermal unbinding transition is called DNA denaturation (see 2] and references therein). It has been modeled 
in various ways, the most prominent being the Poland-Scheraga (PS) model 2lj. Even though this model does 
not take into account spatial aspects of the denaturation transition, it correctly treats sequence effects, and has 
been numerically implemented in the program MELTSIM 0. It has also been shown that a mechanically induced 
denaturation transition is possible, by the combined application of a stretching force and an untwisting torque (see 
and references therein). In a brief review of the standard PS model, we emphasize the role of the loop exponent 
c on the existence and nature of the denaturation transition. We generalize the theory to include a torque which 
introduces a torsional enthalpy term, and results in a modified loop exponent c' < c. Numerical simulations on 
biological sequences show that the denaturation curves are shifted, while the peaks are smoother than their zero torque 
counterparts. Extension to supercoiled or undercoiled DNA , through the use of a numerical Legendre transform, 
yields denaturation isotherms for the same sequences. 

We briefly review the Poland-Scheraga (PS) model for DNA melting, and consider a double stranded (ds) DNA 
fragment, made of N complementary base pairs, assuming that bases (1) and (N) on both strands are paired. We 
denote by Z{a) the forward partition function of the two strands, starting at base (1) and ending at base (a), with 
bases (a) being paired. This partition function satisfies the recursion relation (Figure 1) 
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where P = l/fc^T is the inverse temperature, Sa.a+i is the stacking energy of base pairs {a, a + 1), and as is the bare 
loop formation (cooperativity) parameter (we assume that as is base independent). 
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FIG. 1: Recursion relation for Z{a + 1) (eq.Q) in the PS model. 

The factor J\f{a'; a + 1) counts the number of conformations of a pair of chains starting at base pair (a') and ending 
at base pair (a + 1). It also represents the number of conformations of a closed polymer of 2(a — a') monomers, which 
is asymptotically given by 
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g{a~ a') 



(a — a'Y 



(2) 



where fc^log/ig is the entropy per base pair (assumed to be independent of the chemical nature of the pair), and 
g{x) = ^ is the probability of return to the origin of a loop of length 2x. The exponent c depends on the interaction of 
the loop with itself and with the rest of the chain: It has been extensively discussed in the context of homopolymcric 
DNA [^jHQ- If O'^'-' neglects the interaction with the rest of the chain, we have c — iv (yielding 3/2 for a Gaussian 
loop and « 1.8 for a self avoiding loop). Taking into account the interaction with the rest of the chain is a difficult 
problem: approximations and numerical calculations point toward a value « 2.15 for the full problem |9lll0l|. 

The recursion relation ^ is supplemented by the boundary conditions Z{\) = 1; Z{2) = e~^^^-'^ Z{1). This 
recursion relation can easily be solved analytically if one assumes that all stacking energies are equal. One may for 
instance introduce a grand canonical partition function Z{z) = Z{a) . We summarize the results of this 

homopolymcric study: i) If 2 < c, there is a first order (discontinuous) unbinding transition, ii) If 1 < c < 2 , there 
is a second order (continuous) unbinding transition, with a specific heat exponent a — ^r5f ■ iii) If c < 1 , the two 
strands are always bound and loops open in a continuous way. 

For non homogeneous sequences, the calculation cannot be done analytically. However, the results pertaining to 
the existence of an unbinding transition are expected to hold. In addition, in order to calculate the probability of 
opening of a base pair, it is necessary to introduce forward and backward partition functions The forward 

partition function Zf{a) is nothing but Z[a) whereas the backward partition function Zb{a) is the partition function 
of the two strands, starting at base {N) and ending at base (a), with base (a) being paired. These points will be 
discussed in detail in [l^ . 

We now generalize the Poland-Scheraga model to the case where a torque is applied to the DNA fragment. 

We again assume that base pairs (1) and (N) of the DNA fragment are kept fixed and apply a weak torque F on 
base pair (N). By "weak", we mean that there are no plectonemes on the chain: Experimentally, applying a force 
F > 0.5 pN on a DNA fragment of a few persistence lengths Ip {Ip « 150 bp in (ds) DNA) is enough to prevent the 
formation of plectonemes [5j . 
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FIG. 2: Schematic representation of the force F + torque F experiment. Basepair (ot) has orientation (Qa). The natural twist 
of the hehx is (So)- 

If the force F is applied along the z— axis, the DNA fragment will be aligned (on the average) along this direction. 
One may then assign, in the xy plane, an angle Qa to a paired base pair (a), representing the angle of this pair with 
the X— axis (fig. 2). We denote by the natural twist angle of the DNA helix per base pair (6*0 = 27r/10.4 in radians), 
and model the torsional energy between neighboring base pairs by 

eo(a,a + 1) = ^Ko(6la+i - 6*^ - 6*0)^ (3) 
where kq is the elastic torsion constant of (ds) DNA. 
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Generalizing the PS model, one may define a partition function Z{a, 9a) of the two strands, starting at base pair 
(1) with orientation 0i, and ending at base pair (a), with orientation Oa- This partition function satisfies 



where e^.a+i again denotes the stacking energy of base pairs (a, a + 1) and {Mq, Mi) are normalization factors (see 
below) . 

In equation the existence of the stretching force F is implicit. The functions Ei(^a' ,9a'',a + l,9a+i) and 
Af{a' , 9a' ; a + 1, 9a+i) represent respectively the torsional energy and the number of conformations of a pair of chains 
starting at base pair (a') with orientation {9a>) and ending at base pair {a + 1) with orientation {9a+i)- 

We first discuss M{^a' , 9a'\a-\- 1, 6*0+1)- With unconstrained orientations, equation Q emphasizes the importance 
of the loop exponent c, representing the interaction of the loop with itself and with the rest of the chain. In a 
forthcoming paper , we show that c can also be modeled through the introduction of a specific repulsive potential 
in an otherwise Gaussian loop. In a nutshell, the two strands partition function can be factorized as a product of a 
center of mass partition function Zcm{F) and a relative coordinate partition function 



Zp = / Vp{s) e 



-I3H, 

with (3Hp = ^ J ds p (s) + D J ds ^^^^|;^ , where a « 50 A is the single strand (ss) DNA Kuhn length. The effective 
chain described by p{s) has length 2 (a — a'), and the probability of return to the origin is given by 

~ - Jdr<r\e-^ic-^'W>^^)\0> 

with phD = -^'^l + -S,. 

This model will mimick the original problem, if one has goic^ — a') ^ (a-^a'yo ■ ^ scaling argument then implies 
X — 2. The precise relation between the strength D of the repulsive potential and the loop exponent cd depends 
on the behaviour of the density of states of the Hamiltonian (3hD at low energy. It can be obtained numerically ; in 
particular D can be chosen so that one has cd = c (e.g. c = 1.8 is recovered with ^ ~ i). 

On the other hand, the orientation constraint can be written as 

A9 = 9a+i - 9a' = (6) 
J a' + y' 

where (a;, y) a re the coordinates p perpendicular to the stretching force F. This type of constraint arises in entangled 
polymers |14j . 

Following a calculation of Wiegel , the use of a directed approximation for the stretched strands enables us to 
write the number of conformations A^(a', 9a'': o + 1, ^a+i) as (compare with eq. (0)) 

AA(a', 9a';a + l,9a+i) - ^^^^ ° hi9a+i - 9a') (7) 

0^ 

where p{F) = Pq e 12 and h{9 a+i — 9a') is a (normalized) measure of the torsional entropy reduction given by 

with A = .^/^^log(|^(a — a')). In the previous equation, d denotes the diameter of the double helix (d k, 20 A) 
and as seen above, a denotes the Kuhn length of (ss) DNA (a w 50 A). The experimental value of log/ip is taken as 



12.5 (see ref.0|). For forces F of order 1 pN, the difference between log/i and log/iQ is of order 0.1 and will thus be 
neglected in the following. The validity of the directed approximation for this calculation will be discussed thoroughly 
in a forthcoming paper |l3| . 

In equation the torsional energy -Ei(a', 9a'; a + I, 9a+i) of a bubble {a' ,a+ 1) is assumed to be of the form 

J7 f ' a I 1 A Ki {9a+l - 9a')'^ , , 

Ei[a ,9a';a + l,9a+i) ^ — — -, — (9) 

2 (a — a ) 

where ki is the torsional constant of a DNA bubble. 

For long enough loops (a — a' >> log(a — a')), this energy is small compared to the entropy reduction of eq. 
Furthermore, due to the softness of unbound fragments, one expects that k; << kq. In the following, we will thus set 
this torsional energy to zero. 

The recursion relation for the partition function therefore reads 

Z{a + l,9a+i)= e-'^--^^ ^e-'-^io.^^-o.-Oof z{a,9a) 



+ tTs^ / -^^{c^',9a';a + l,9a+i)Z{a',9a') (10) 

where A/'(a', 0q'; a + 1, 6*0+1) is given in eqs. (|3|H|). 

Since the integration over the angular variables 9a should yield back equation the normalization factors are 

Mo = ^ and M, = 1. 

Setting 9i = 0, the boundary conditions pertaining to equation (|10() are Z{l,9i) — S{9i) and Z{2,92) — 
j_„-/3£i,2-^(e2-eo)" 

Equation IjlOfl can be brought to the form of a standard PS recursion relation by going to the torque representation. 
We define the Laplace transform 

/+00 
d9ae'^^'^''Z{a,9a) (11) 
-00 

The quantity Z(a, F) represents the partition function of a DNA chain of length a fixed at the origin, and subject 
to a torque F. Taking the Laplace transform of H10() . we obtain 



Z{a + 1, F) = e-'^'^'-.o+i Z{a, T) + a'g 



Mo 



a'=i (a - 



- Z(a',F) 



(12) 



with 



r9n 



2ko 



(13) 
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The boundary conditions translate into Z{1,T) = 1 and Z{2,T) = e'^"^^-^ 

In the form (|12|l . we recognize a standard PS recursion equation, with stacking energies given by H13(l . loop exponent 
given by H14|) and loop formation (cooperativity) parameter given by (|15|l . These new effective parameters have the 
following properties: i) the loop exponent is decreased by the torque, so that the probability of return to the origin 
is increased (see eq.Q), ii) the loop formation parameter is increased by the torque. With realistic values of the 
parameters (see below), this effect is very weak. 

These results can easily be understood for undercoiling (F < 0). However, in the case of strong positive supercoiling 
(r > 0), the spatial arrangements of the bases are important, a feature which is absent of the PS approach. This can 
lead to the appearance of new phases such as P-DNA j^. We thus can trust our approach for negative and weakly 
positive torques. 

If one takes all stacking energies equal to —Eq , equation H12|) can be solved analytically by using the appropriate 
grand canonical partition function Q{z,r) = z"Z{a,T) with the results Q{z,r) = - 



— jp — where 

Pc'{z) = p"- As in the original PS model, the critical point is obtained when z is a pole of the denominator 

and z/iQ = 1. As previously mentioned, the specific heat exponent (in the thermodynamic limit) reads a{T) — ^§3p. 

A theoretical outcome of our model is that the denaturation transition disappears at large enough torque. Indeed, 
as we perviously saw, the PS model displays a phase transition only if the loop exponent c is larger than 1. Equation 
()14|l shows that even though the bare exponent c is larger than 1, it becomes smaller than 1 when the torque is 
increased. Therefore, the transition gets smoothed out as the torque is increased: The denaturation peaks broaden 
and are shifted to lower or higher temperatures, depending on the sign of the applied torque. 

For non homogeneous stacking energies, the properties of Q{z,T) are not amenable to analytic calculations and 
one has to resort to numerical calculations. The parameters we use are c ~ 1.8 (corresponding to Dja} k, 0.33), 
OS = 1.26 10-^ and the MELTSIM stacking energies 

In Figure 3, we plot the derivative with respect to the temperature of the fraction of bound pairs c — —^ (related to 
the experimental UV absorption of DNA) for a biological sequence 0| of 2000 base pairs as a function of temperature, 
for various values of the torque F. As F increases, the peaks get smoothed out. 
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FIG. 3: Derivative of the number of paired bases c = — 
{thin), -150 (thick), -300 (thin), -500 (thick) K, from right to left. 



for a 2000 bps biological fragment, for F — 



Up to here, we have considered the problem of a DNA fragment subject to a weak torque. We now come back 
to the case where the DNA fragment is supercoiled or undercoiled (as is the case of circular DNA in plasmids): the 
winding angle 0jv ~ is not equal to the total natural twist angle (N — 1)0q . The supercoiling index s is defined by 



^jv-^i = (N-l) 9o (1 + s) 



(16) 
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It can be positive (supercoiling, F > 0) or negative (undercoiling, F < 0). 

The partition function Z(N,9]\j) is related to the torque representation through an inverse Laplace transform 



/•Co+ioo jp 

Z{N,0n)^P ^e-^^'^^-'^'>°^'+'^ Z{N,T) (17) 

Jco-ioo 2i7r 

where Cq is a constant which leaves all the singularities of Z{N,T) to its right. For large N, one may perform a 
saddle point calculation, with F^ defined by 

1 91ogZ(7V£) 

J^i )r=r, = /30o(l + s) (18) 



For homogeneous sequences (all stacking energies equal to — £o), one may go one step further since 



ZiN,0^) ^pi r^^^ ^^-pnN-,)e.ii+s) Q(^^r) (19) 



where the z integral is to be performed on a circle containing the point 0. 

In the limit of large iV, equation H19() can again be evaluated by the saddle point method on both z and F . The 
saddle-point equation for z is given by = z^gf-logQ . In the thermodynamic limit, N —> oo, the saddle-point 
solution Zs should approach the (F-dependent) pole z* of Q{z,T), defined by 



The saddle point H18|) then reads 



and to leading order, the free energy per base pair /(s) is given by 



1 - z*e'^'''° - (t',z* Pc' {z*fi„) = (20) 



/3^o(l + .) = (^W (21) 



/3/(s) = logz*(F,)-f/3F,0o(l + s) (22) 

To make contact with previous work, we restrict ourselves to small supercoiling (s << 1) or equivalently small 
couple F. U sing equation (|20|l we can expand z* to second order in F . The details of the calculations will be 
presented in '13] and we just quote the results: 



(3m = log.*(F = 0) + ^ + (23) 

where A2 is a constant, equal to the fluctuation of the end angle (A2 = gr^^ )r-o ~ ^{(^n) ~ (^^')^)r-o) 

The physical interpretation of equation (|23|l is fairly simple. It essentially states that for small supercoiling s, the 
free energy of a DNA strand is equal to the sum of the free energy of the non supercoiled fragment plus a free energy 
term which forces the fraction of unbound pairs to be close to the opposite of the supercoiling index. This form is very 
similar to the form which was devised phenomenologically by Benham In our case, this form is derived from the 
microscopic model, and the coefficients of the quadratic part are expressed in terms of the microscopic characteristics 
of the DNA. 

For larger supercoiling, the relation between s and the fraction of unbound pairs is not straightforward. Arguments 
that suggest that — s ~ for any s are given in ref. 

Finally, we point out that the specific heats C(T, F) and C{T,s) are related through Fisher renormalization [2(tI |. 
In particular, if C(T, F) diverges, C(T, s) will be finite at the transition. 

For non homogeneous sequences, numerical calculations yield Z{N, F), and through ea. H18|) . one obtains the isotherm 
curves F(s). The application to the sequence used in Figure|31is shown in Figure^] 
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FIG. 4; Torque as a function of the supercoiling index for various temperatures, for the sequence of Figure 3. 



The statistical mechanics of the torque induced DNA denaturation has been previously studied in different contexts 
[MI^HMH^i and there is agreement with our results whenever they overlap. Ref 21] considers a transfer matrix 
formalism for homogeneous sequence with c = |. Ref starts from the Benham phenomenological expression and 
obtains - among other results - the critical curve Tc(s) in the homogeneous case. Ref |23| studies the thermodynamic 
and dynamic properties of a random DNA fragment submitted to a torque. Finally Ref |5( discusses the experimental 
procedure in detail (ranges of force and torque); the comparison with our results is not straightforward as experiments 
rely on the extension curve, which is not included in PS type of models. 

We have proposed a generalization of the PS model to include the effect of torsion on DNA denaturation. Torsion 
produces a very strong reduction of entropy in the loops, which eventually suppresses the denaturation transition. 
Application to cyclic DNA (plasmids) will be the subject of future work. 
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